Momentum Distribution and Contact of the Unitary Fermi gas 
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We calculate the momentum distribution n{k) of the Unitary Fermi Gas using Quantum Monte 
Carlo calculations at finite temperature T/ep as well as in the ground state. At large momenta k/kp, 
we find that n{k) falls off' as C/k'^, in agreement with the Tan relations. From the asymptotics of 
n(fc), we determine the contact C as a function of T/cp and present a comparison with theory. 
At low T/cp, we find that C increases with temperature, and we tentatively identify a maximum 
around T/tp ~ 0.4. Our calculations are performed on lattices of spatial extent up to A^^, = 14 with 
a particle number per unit volume of ~ 0.03 — 0.07. 

PACS numbers: 67.85. De, 67.10.Jn, 05.30. Fk 



I 

o 



(N 
> 

in 



X 



The Unitary Fermi Gas (UFG) is one of the most in- 
teresting strongly interacting systems known to date, as 
it saturates the unitarity bound on the quantum me- 
chanical scattering cross section dg < ATr/k"^. Since 
the proposal of the UFG as a model for dilute neutron 
matter by Bertsch [l[ and its realization in ultracold 
atom experiments [1] , the UFG has garnered widespread 
attention across multiple disciplines, including atomic 
physics 3, nuclear structure Q and relativistic heavy- 
ion collisions Q . The UFG is defined as a two-component 
many-fermion system in the limit of short interaction 
range and large s-wave scattering length a, 



<— kpT^ ^ 1 ^ kpa — ?> oo. 
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with kp = (Stt^ti)^/'^ the Fermi momentum and n the 
particle number density. The special properties of the 
UFG arise from the fact that it is characterized by a sin- 
gle scale, given by the inter- particle distance kp^ , with- 
out reference to the details of the interaction. While the 
thermodynamic properties of the UFG are universal Q , 
the lack of an obvious dimensionless expansion parameter 
makes the UFG a challenging many-body problem. 

In spite of the challenges of the unitary limit, much 
progress has been made with purely analytical methods. 
Notably, in 2005 Tan was able to derive exact thermody- 
namic relations Q in terms of a universal quantity known 
as the "contact" C, which determines the number of pairs 
separated by short distances. Since then, the Tan rela- 
tions have been re-derived in multiple ways 0413 as well 
as verified experimentally 
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Recently, C has also been found to determine the pref- 
actor of the high-frequency power-law decay of correla- 
tors [3, 113 ) as well as the right-hand sides of the shear- 
and bulk viscosity sum rules [13 • The contact is there- 
fore a central piece of information on the UFG in equilib- 
rium as well as away from equilibrium, since it constrains 
several thermodynamic quantities with a single number. 



On the experimental side, C has been shown to be cen- 
tral to radio-frequency spectroscopy and laser photoas- 
sociation [l^, as well as to govern the rate of decrease 
of low-energy atoms due to inelastic two-body scattering 
processes with a large energy release. The Tan relations 
(as well as the above-mentioned sum rules) remain valid 
at arbitrary kpa as long as kpT^ <^ 1. For further details 
and a comprehensive review, see Ref. jl6j. 

The calculation of C itself, however, remains a chal- 
lenge, as it depends on the intricate many-body dynamics 
of the unitary regime. In principle, C can be extracted 
from any one of the Tan relations (as recently done in 
experiments [13] ) ■ One of the simplest relations concerns 
the asymptotics of the momentum distribution, and as- 
serts that 



C = lim k*n^{k), 
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where n^{k) is the momentum distribution expressed as 
a thermal average, and the ^ and a^ j. denote creation 
and annihilation operators for particles of momentum k 
and spin a. If n^{k) is normalized to the particle number 
N^, then C is an extensive quantity with dimensions of 
momentum. We shall consider C in units of kp divided 
by the total particle number N — -\- N^. 

In this work, we focus on the momentum distribution 
of the homogeneous UFG and the extraction of C via 
Eq. using a Quantum Monte Carlo (QMC) approach 
which accounts fully for quantum and thermal fluctua- 
tions. On a spatial lattice, the Hamiltonian that captures 
the physics of the unitary limit can be written as 
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where m is the mass of the fermions (henceforth set to 
unity), g is the bare coupling, and n^^ denotes the num- 
ber density operator for spin projection a at lattice posi- 
tion i. The equilibrium thermodynamical properties are 
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obtained from the grand canonical partition function 

Z = Tr exp[-/3(i/-^iV)], (4) 
where (3 = l/k^T, fi is the chemical potential, and 



Hi 



(5) 



denotes the particle number operator. 

In our QMC treatment, the system is placed on a 
(3 + l)-dimensional Euclidean space-time lattice via a 
Suzuki- Trotter decomposition of the Boltzmann weight 
in Eq. ([4]), and the interaction is represented via a 
Hubbard-Stratonovich (HS) transformation [l7|. As we 
focus on the spin-symmetric case, the fermion sign prob- 
lem is absent. The resulting path integral formulation 
is an exact representation of the many-body problem of 
Eq. (Ul, up to finite volume and discretization effects. 
These may be addressed by varying the spatial lattice 
volume V = iV^ and the density n, such that the thermo- 
dynamic and continuum limits are recovered as — ;> oo 
and n 0. The latter requires great care, as too low den- 
sities imply a departure from the thermodynamic limit. 
We find that n ~ 0.03 — 0.05 particles per unit volume 
yield results accurate to ~ 7% at finite temperature, and 
to < 5% at T = 0. 

Our lattice formulation is very similar to Ref. jl^ j but 
differs in at least three notable aspects. Firstly, we deter- 
mine the bare lattice coupling constant g corresponding 
to the unitary regime by using Liischer's formula [l9j as 
in Ref. (l^. This procedure yields g ~ 5.14 in the uni- 
tary limit. Secondly, we use the compact, continuous HS 
transformation 



exp [rg n^,n^.^ 
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where (not to be confused with the spin projection) is 
the HS auxiliary field, with /2 = cxp(T5) — 1, and t de- 
notes the lattice spacing in the imaginary time direction. 
We find that r ~ 0.05 is sufficiently small to render dis- 
cretization errors from the Suzuki- Trotter decomposition 
insignificant (see also Fig. [5]). The above representation 
(referred to as "Type 4" in Ref. [2l|) was found to be su- 
perior with respect to acceptance rate, decorrelation and 
signal-to-noise properties than the more conventional un- 
bounded and discrete forms ji^l . Thirdly, we update the 
auxiliary field a by using the Hybrid Monte Carlo (HMC) 
algorithm [23)] (familiar from Lattice QCD), which com- 
bines the Metropolis algorithm with deterministic Molec- 
ular Dynamics. Our implementation of the HMC algo- 
rithm enables global updates at all temperatures and lat- 
tice sizes, and scales approximately as ~ as a function 
of the spatial lattice volume, to be contrasted with the 
~ scaling of approaches based on local updates. 




FIG. 1: (Color online) Momentum distribution n{k) from 
QMC for A*'^ = 10 as a function of k/kp, for various tem- 
peratures ranging from zero to T/ep ~ 0.5. The curves are 
intended as a guide to the eye, and the statistical errors are 
the size of the symbols. Inset: n{k) for = 14 in a log-log 
scale, showing the asymptotic ~ behavior. 



We have performed calculations at T = as well as 
T/ep > 0, in the former case using an approach similar 



to Ref. [21| . Our main results correspond to 40 — 50 par- 
ticles at = 10 and 70 - 80 particles at = 12, in 
addition to limited data for — 14. In Fig. [1] we show 
the momentum distribution n{k) as a function of T/ep. 
We have computed n(k) by averaging over the angular 
directions on the lattice as well as over the imaginary- 
time slices. In this way, we find that ^ 200 uncorrelated 
auxiliary field samples for each datapoint gives excellent 
statistics for n{k). Multiplying n{k) by k'^, as plotted 
in Fig. [21 we find a peak at k ^ kp and a leveling out 
at high momenta, with the asymptotic regime setting in 
at fc ~ 2kp at the lowest temperatures. It is fortuitious 
that the asymptotic regime sets in at such low momenta, 
as there is no obvious reason for this to be the case. It 
is then possible to study the temperature dependence of 
this "plateau", which allows for a determination of the 
contact C/{Nkp) as a function of T/ep. These results 
are given in Fig. [31 together with a comparison with avail- 
able theoretical analyses. Our results indicate that n{k) 
follows the expected ^ /c"^ dependence very accurately 
up to at least k ~ 4fcj^, at which point the signal deteri- 
orates due to lattice artifacts. 

The value of C in the ground state can be computed 
via Diffusion Monte Carlo (DMC) calculations, as first 
done in Ref. [i^l using density-density correlations, which 
yielded C(T = 0)/{Nkp) ~ 3.4, up to errors associated 
with fixing the nodes of the wavefunction. A more re- 
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N„ = 1 2, T/Ef = 0.404, p = 0.033 
N, = 12, T/Ep = 0.178, p = 0.040 
C/Nk, = 3.32 
C/Nk, = 3.08 



2 2.5 

k/kp 



4.5 



4.5 - 



3.5 



Z 

o 



2.5 



T-matrix, Palestini et al. 
T-matrix, Enss et al. 




0.2 0.4 



0.6 

T/er 



1.2 



FIG. 2: (Color online) Plot of ?,TT'^{k/kpfn{k) for = 12 
as a function of k/kp at T/ep = 0.178 and 0.404. The 
"plateaux" at large k/kp give the (intensive) dimensionless 
quantity C/(Nkp). At low T/tp, the asymptotic region is 
reached at k/kp ~ 2. Inset: A^^ = 10 results at T = show- 
ing only slight dependence on the Suzuki- Trotter step r. 



FIG. 3: (Color online) Summary of QMC resuhs for C/{Nkp) 
as a function of T/ep, as determined from the large k/kp be- 
havior of n{k). The errorbars are dominated by systematics 
related to the residual fluctuations in the plateaux, as exhib- 
ited in Fig. [21 Also shown are the t-matrix calculations of 
Ref. [27! . [2^ . the virial expansion of Ref. [2^ and the dia- 
grammatic Monte Carlo result of Ref. [33 |. 



cent and comprehensive DMC calculation [25| came to 
the same conclusion using the equation of state, the 
momentum distribution and the density-density corre- 
lation. In contrast, our present results indicate that 
C{T = 0)/{Nkp) ~ 2.95 ± 0.10. The cause of this dis- 
agreement is being explored. The main sources of uncer- 
tainty in our determination of C/{Nkp) are due to finite 
density effects. While we find that such effects tend to 
overestimate C/{Nkp) as well as degrade the formation 
of an asymptotic ^ tail in n{k) at larger values of 
T/ep, larger lattices are needed in order to maintain the 
thermodynamic limit at lower densities. 

The temperature dependence of C at unitarity was 
first determined analytically in Ref. [26|, who consid- 
ered two different limits. At very low temperatures 
T <C T|, ~ 0.15 e^, the dominant excitations are of 
phononic origin, and the T-dependence of C is of the form 
C/{Nkp) oc {T/ep)'^. On the other hand, at very high 
temperatures T:$>ep, one finds C/(7Vfcj,) ~ 16/3 (e^/T) 
within the second-order virial expansion. An interpo- 
lation between these limits then suggests that C{T/ep) 
should present a maximum for T ^ Cp. Recently, C has 
also been computed using two different types of t-matrix 

as well as a third-order virial 



approximations 
expansion (29j . 
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The latter has shown evidence for con- 



vergence of the virial expansion down to T ^ e^. In light 
of these findings and upon analysis of various model cal- 
culations at low T, Ref. [2^ conjectured that the contact 
is likely a monotonically decreasing function of T, except 
possibly in the phononic regime at very low T . While 



the virial expansion is on solid ground at high T, where 
it agrees with the t-matrix approaches of Refs. 
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the actual T-dependence in the strongly correlated low-T 
regime has remained an open question, particularly since 
the UFG is strongly correlated even above [sOl ■ 

Our results show that C grows with T well beyond 
the superfluid phase, and are suggestive of a maximum 
C^jjx ~ 3.4 at T/ep ~ 0.4. This scenario is in quali- 
tative agreement with Ref. (2^, as well as the t-matrix 
calculation of Ref. [13] ■ As C measures the number of 
particle pairs (of both spins) whose separation is small, 
the appearance of a maximum indicates an enhancement 
in such short-range correlations. This may be a result 
of local pairing order which in turn suggests that 
C'max is directly related to pairing above T^, i.e. to a 
pseudogap. We find the scale at which the k~'^ law sets 
in (see Fig. [2]) to be /c ~ 2kp at finite T/ep and somewhat 
lower for the ground state, in agreement with Ref. 12 1. 
This universal property of the unitary limit character- 
izes the "healing distance" of the two-particle boundary 
condition on the many-body wavefunction, and therefore 
separates the microscopic properties from the universal 
macroscopic aspects of the unitary regime. Direct com- 
parison of our data with ultracold atom experiments can 
be achieved by means of the virial expansion and the Lo- 
cal Density Approximation (LDA). While we defer this 
issue to a follow-up paper, wc note that in light of the 
work of Ref. [31}, the features of C(T/ep) found in this 
study are unlikely to conflict with current experiments. 
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In summary, we have computed the momentum dis- 
tribution n{k) and the contact C/{Nkp) for the UFG 
at zero and finite T/ep, using the auxihary field QMC 
method in conjunction with the HMC algorithm. While 
the ground-state momentum distribution was first de- 
termined via DMC calculations in Ref. [32], our results 
represent the first fully non-perturbative calculation of 
n{k) free of uncontrolled approximations. We find that 
the contact at T = assumes the value ~ 2.95 ± 0.10 
and increases as a function of T/tp in the low- and in- 
termediate temperature regimes that we have explored, 
which is consistent with the phononic scenario. No- 
tably, DMC calculations find a somewhat larger value of 
C/{Nkp) ~ 3.4, while the analytic approach of Ref. [33| . 
which interpolates smoothly between the strong- and 
weak-coupling limits, yields C/{Nkp) ~ 3.0 which is con- 
sistent with our data. Our results complement the calcu- 
lations of Refs. (26l - [29| . and are suggestive of a maximum 
in C/(Nkp) at T/ep ~ 0.4, which agrees qualitatively 
with Ref. ^ but disagrees with Ref. [H]. While cal- 
culations at higher T/ep ^ 1 are feasible, an improved 
understanding of the finite density effects is clearly called 
for. 
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